---
title: "McCourt and Ruley Analysis Replication"
author: "Garrett Ruley"
date: '2022-06-10'
output: word_document
---

# Readme

This document contains the code used to create the analysis and visualizations used by McCourt and Ruley in the paper "How is the Establishment Structured? A Multiple Correspondence Analysis of the U.S. China Field."  The code has been broken into chunks for each stage of the analysis, and each chunk has short explanation in plain text at the start as well as comments throughout.

To replicate our results, the only step the reader must take is to point the first command in the third chunk - a command to load the data from the local drive into R - to the correct file.  Running all chunks will then produce our results and their visualizations.  Figures and tables will be stored in the R environment and can be written to the local drive if desired.

# Replication Code

Load necessary packages
```{r setup and packages, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(readxl) #For loading Excel data files
library(tidyverse) #For data cleaning and manipulation
library(knitr) #For knitting markdown to Word
library(xtable) #For creating clean tables of results
library(FactoMineR) #For conducting MCA
library(ggplot2) #For visual plotting of MCA results
library(ggrepel) #Additional plotting tool for embedding text in figures.
library(sqldf) #Allows use of SQL commands for data management
```

Define custom convenience functions that aid in creating tables of key statistics.
```{r EZMCA, include = FALSE}
#This function calculates "modified rates of variance" for each of the principal axes in an MCA.  This method of calculating modified rates is reproduced from Le Roux and Rouanet (2010), 

modified.rates <- function(MCA) {
 
  Q <- length(names(MCA$call$X))
  Q <- as.numeric(Q) 
  
  lambbar <- 1/6
  
  R <- as_tibble(MCA$eig) 
  
  lmax <- subset(R, eigenvalue>lambbar) 
  lmax <- lmax %>% pull(eigenvalue) 
  
  i <- length(lmax)
  i <- i[1]
  i <- as.numeric(i)
  
  pseudoeig <- numeric(i) 
  
  for(val in c(1:i)){
    pseudoeig[val] <- ((6/(6-1))^2)*((lmax[val]-lambbar)^2)
  } 

  S <- sum(pseudoeig) 
  
  modrates <- pseudoeig/S 
  modrates
}

#This function produces a basic table of the modified rates, and the cumulative modified rates, by axis.

mod.table <- function(x){
  modrates <- x
  
  l <- length(x)
  
  summod <- numeric(l)
  
  for(val in 1:l){
    partsum <- modrates[1:val]
    
    summod[val] = sum(partsum)
  }
  
  modtable <- cbind(modrates, summod)
}
```

Load dataset and conduct MCA.  *This is the chunk where the reader must ensure that the "read.csv" command points to the correct address for our data on the reader's machine.*
```{r run MCA, include = FALSE}
basedata <- read.csv("INSERT FILE PATH HERE")
basedata$X <- NULL

#The "shortname" variable will be useful in labeling graphs later, but clearly cannot be used as either an active or a passive variable in the MCA.
MCAdata <- sqldf('SELECT orgtype, catdistance, catfirstyear, partisanalignment, catnumspecialists, catscore, catloc, catobama, cattrump FROM basedata')

#After all that, the simple FactoMineR command MCA() runs the analysis.  By specifying "quali.sup = 6,8,9", the sixth, eigth, and ninth variables in the dataset ("catscore", "catobama", and "cattrump") are coded as passive.
ROCMCA <- MCA(MCAdata, ncp = 4, quali.sup = c(6,8,9))
```

Prepare tables.  These tables contain both base statistics and the MCA results used to create visualizations.

Tables with data included directly in the paper:
basestats --> 
baseorgtable --> 
basepartable --> 
basedistable --> 
baseloctable --> 
basespectable --> 
baseyeartable --> 

Tables with data used for visualization:
tidyvartable --> 
tidyindtable --> 
tidyvar1 --> 
tidyvar2 --> 
tidyvar3 --> 
tidyind1 --> 
tidyind2 --> 
tidyind3 --> 
twittertable --> 
wholetwittertable --> 
speechtable --> 
wholespeechtable --> 
```{r table prep, include = FALSE}
#Create tables of coordinates and contributions for both the cloud of categories and the cloud of individuals.  Data is regularly extracted from the ROCMCA object, which contains the raw results of the MCA but in a difficult-to-graph format.
varnames <- dimnames(ROCMCA$var$coord)
varnames <- varnames[1]
varnames <- unlist(varnames)
  
tidyvarcoord <- as_tibble(ROCMCA$var$coord)
names(tidyvarcoord)[1:4]<-c("Axis1", "Axis2", "Axis3", "Axis4")
  
tidyvarcontrib <- as_tibble(ROCMCA$var$contrib)
names(tidyvarcontrib)[1:4] <- c("Contrib1", "Contrib2", "Contrib3", "Contrib4")
  
tidyvartable <- as_tibble(cbind(varnames, tidyvarcoord, tidyvarcontrib))
  
#Individuals

tidyindcoord <- as_tibble(ROCMCA$ind$coord)
names(tidyindcoord)[1:4] <- c("Axis1", "Axis2", "Axis3", "Axis4")
  
  tidyindcontrib <- as_tibble(ROCMCA$ind$contrib)
  names(tidyindcontrib)[1:4] <- c("Contrib1", "Contrib2", "Contrib3", "Contrib4")
  
  tidyindtable <- as_tibble(cbind(tidyindcoord, tidyindcontrib))

indnames <- basedata %>% pull(shortname)
tidyindtable <- as_tibble(cbind(indnames, tidyindtable))

#Pull the coordinates of the supplementary variable into a separate table
steptable <- tidyvartable[,1:5]

suptable <- ROCMCA$quali.sup$coord
supnames <- dimnames(suptable)
supnames <- supnames[1]
varnames <- unlist(supnames)

supcoord <- as_tibble(ROCMCA$quali.sup$coord)
names(supcoord)[1:4]<-c("Axis1", "Axis2", "Axis3", "Axis4")

suptable <-  cbind(varnames, supcoord)
wholetable <- rbind(steptable, suptable)

#Base response patterns

#Set our important constants
K <- 28 #Number of categories
N <- 87 #Number of responses
Q <- 6 #Number of questions

#Create a table of base statistics (number, relative frequency, and category contribution).  Formulas taken from Le Roux and Rouanet (2010).
basestats <- data.frame(n = numeric(28),
                      f = numeric(28),
                      Ctr = numeric(28),
                      stringsAsFactors=FALSE)

# "For" loops are used to compute relative frequencies for the categories within each variable
for(val in 1:28){
  basestats$n[val] <- sum(ROCMCA$call$Xtot[,val])
  basestats$f[val] <- basestats$n[val]/N
  basestats$Ctr[val] <- ((1-basestats$f[val])/(K-Q))
}

varnames <- dimnames(ROCMCA$var$coord)
varnames <- varnames[1]
varnames <- unlist(varnames)

basestats[,1:3] <- round(basestats[,1:3], digits = 3)

basestats <- basestats %>% mutate(f = f*100)
basestats <- basestats %>% mutate(Ctr = Ctr*100)

rownames(basestats) <- varnames
colnames(basestats) <- c("n", "Frequency (in %)", "Contribution")

#Create question-specific tables
#Orgtype

baseorgtable <- basestats[1:8,]
baseorgtable[9,] <- c(0,0,0)

baseorgtable[9,1] <- sum(baseorgtable[,1])
baseorgtable[9,2] <- sum(baseorgtable[,2])
baseorgtable[9,3] <- sum(baseorgtable[,3])
rownames(baseorgtable)[9] <- c("**Totals**")

#Partisan Alignment

basepartable <- basestats[20:22,]
basepartable[4,] <- c(0,0,0)

basepartable[4,1] <- sum(basepartable[,1])
basepartable[4,2] <- sum(basepartable[,2])
basepartable[4,3] <- sum(basepartable[,3])
rownames(basepartable)[4] <- c("**Totals**")

#First Year

baseyeartable <- basestats[14:19,]
baseyeartable[7,] <- c(0,0,0)

baseyeartable[7,1] <- sum(baseyeartable[,1])
baseyeartable[7,2] <- sum(baseyeartable[,2])
baseyeartable[7,3] <- sum(baseyeartable[,3])
rownames(baseyeartable)[7] <- c("**Totals**")

#Number of specialists

basespectable <- basestats[23:26,]
basespectable[5,] <- c(0,0,0)

basespectable[5,1] <- sum(basespectable[,1])
basespectable[5,2] <- sum(basespectable[,2])
basespectable[5,3] <- sum(basespectable[,3])
rownames(basespectable)[5] <- c("**Totals**")

#Distance from capitol

basedistable <- basestats[9:13,]
basedistable[6,] <- c(0,0,0)

basedistable[6,1] <- sum(basedistable[,1])
basedistable[6,2] <- sum(basedistable[,2])
basedistable[6,3] <- sum(basedistable[,3])
rownames(basedistable)[6] <- c("**Totals**")

#Number of US locations

baseloctable <- basestats[27:28,]
baseloctable[3,] <- c(0,0,0)

baseloctable[3,1] <- sum(baseloctable[,1])
baseloctable[3,2] <- sum(baseloctable[,2])
baseloctable[3,3] <- sum(baseloctable[,3])
rownames(baseloctable)[3] <- c("**Totals**")

#Table of variances, variance rates, and modified rates, for principal axes.

axistable <- ROCMCA$eig[1:5,1:2]

modifiedrates <- modified.rates(ROCMCA)
modtable <- mod.table(modifiedrates[1:5])

axistable <- cbind(axistable,modtable)

rownames(axistable) <- c("Dim 1", "Dim 2", "Dim 3", "Dim 4", "Dim 5")
colnames(axistable) <- c("Axis Eigenvalue", "Percentage of Variance", "Modified Variance", "Cumulative Modified Variance")

axistable[,1:4] <- round(axistable[,1:4], digits =4)

#Coordinate and contribution tables
vartable <- tidyvartable
names(vartable) <- c("Question Category", "Axis 1 Coordinate", "Axis 2 Coordinate", "Axis 3 Coordinate", "Axis 4 Coordinate", "Contribution to Axis 1", "Contribution to Axis 2", "Contribution to Axis 3", "Contribution to Axis 4")

#Pull highest-contributing categories and individuals into separate tables for plotting.
tidyvar1 <- sqldf('SELECT * FROM tidyvartable WHERE Contrib1 > (SELECT avg(Contrib1) FROM tidyvartable)') 
tidyvar2 <- sqldf('SELECT * FROM tidyvartable WHERE Contrib2 > (SELECT avg(Contrib2) FROM tidyvartable)')
tidyvar3 <- sqldf('SELECT * FROM tidyvartable WHERE Contrib3 > (SELECT avg(Contrib3) FROM tidyvartable)')
tidyvar4 <- sqldf('SELECT * FROM tidyvartable WHERE Contrib4 > (SELECT avg(Contrib4) FROM tidyvartable)')

tidyind1 <- sqldf('SELECT * FROM tidyindtable WHERE Contrib1 > (SELECT avg(Contrib1) FROM tidyindtable)')
tidyind2 <- sqldf('SELECT * FROM tidyindtable WHERE Contrib2 > (SELECT avg(Contrib2) FROM tidyindtable)')
tidyind3 <- sqldf('SELECT * FROM tidyindtable WHERE Contrib3 > (SELECT avg(Contrib3) FROM tidyindtable)')
tidyind4 <- sqldf('SELECT * FROM tidyindtable WHERE Contrib4 > (SELECT avg(Contrib4) FROM tidyindtable)')

#Separate supplemental variables for plotting
twittertable <- suptable[1:4,]
obamatable <- suptable[5:7,]
trumptable <- suptable[8:10,]
speechtable <- suptable[5:10,]
wholespeechtable <- rbind(steptable, speechtable)
wholetwittertable <- rbind(steptable, twittertable)

Administration <- c("Obama","Obama","Obama","Trump","Trump","Trump")
speechtable <- cbind(speechtable, Administration)

```

Create visualizations using ggplot2 and the tabulated data from the previous section.  The first graph has been completely commented out for reference.
```{r}
#Set ticks and limits for code readibility:
ticks = c(-2.0, -1.5, -1.0, -0.5, 0, 0.5, 1, 1.5, 2.0)
lims = c(-1.5, 1.5)
lims2 = c(-1.5, 2.0) # Depending on the axes selected, some plots need slightly longer axes than others

#Graph categories on principal axes 1-2
plot1.1.12 <- ggplot(tidyvartable,
                   mapping = aes(x = Axis1, y = Axis2, label = varnames)) + # Pull data
  geom_vline(xintercept = 0, color = "gray80") + # Add vertical and horizontal axes to plot
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(alpha = 0.7) + # Plot the categories as points on axes 1-2
  geom_text_repel(size = 2) + # Add text labels with the category names
  coord_fixed(ratio = 1) + # Scale the x- and y-axes equally with each other
  scale_x_continuous(breaks = ticks, limits = lims2) + # Specify tickmarks and axis length for x
  scale_y_continuous(breaks = ticks, limits = lims) + # Specify tickmarks and axis length for y
  labs(x = "Dimension 1 (10.23%)", #Add labels and title
       y = "Dimension 2 (9.80%)",
       title = "Cloud of Categories",
       subtitle = "Principal Axes 1 and 2")

#Graph categories on principal axes 2-3
plot1.1.23 <- ggplot(tidyvartable,
                   mapping = aes(x = Axis2, y = Axis3, label = varnames)) +
  geom_vline(xintercept = 0, color = "gray80") +
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(alpha = 0.7) +
  geom_text_repel(size = 2) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Cloud of Categories",
       subtitle = "Principal Axes 2 and 3")


#Graph the cloud of individuals on principal axes 1-2
lims3 = c(-1.25, 1.80)

plot1.2.12 <- ggplot(tidyindtable,
                   mapping = aes(x = Axis1, y = Axis2, label = indnames)) +
  geom_vline(xintercept = 0, color = "gray80") +
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(alpha = 0.6, size = 1.5) +
  geom_text_repel(size = 1.8) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims3) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  labs(x = "Dimension 1 (10.23%)",
       y = "Dimension 2 (9.80%)",
       title = "Cloud of Individuals",
       subtitle = "Principal Axes 1 and 2")

plot1.2.12
ggsave("COI - Axes 1-2.png", height = 8, width = 8, device = "png")

#Graph the cloud of individuals on principal axes 2-3
plot1.2.23 <- ggplot(tidyindtable,
                   mapping = aes(x = Axis2, y = Axis3, label = indnames)) +
  geom_vline(xintercept = 0, color = "gray80") +
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(alpha = 0.5) +
  geom_text_repel(size = 2) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Cloud of Individuals",
       subtitle = "Principal Axes 2 and 3")

#Categories with contributions above the mean for each principal axis
#Axis 1
plot2.1.1 <- ggplot(tidyvar1,
                mapping = aes(x = Axis1, y = Axis2, label = varnames)) +
  geom_vline(xintercept = 0, alpha = 0.8) +
  geom_hline(yintercept = 0, alpha = 0.8) +
  geom_point(alpha = 0.8, size = 4, mapping = aes(color = Contrib1)) +
  geom_text_repel(size = 3, point.padding = 0.1) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims2) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  scale_colour_distiller(palette = "YlOrRd", direction = 1) +
  labs(x = "Dimension 1 (10.23%)",
       y = "Dimension 2 (9.80%)",
       title = "Main Category Contributors to Axis 1") +
  guides(color = guide_colourbar(title = "Degree of Contribution"))

#Axis 2
plot2.1.2 <- ggplot(tidyvar2,
                mapping = aes(x = Axis2, y = Axis3, label = varnames)) +
  geom_vline(xintercept = 0, alpha = 0.8) +
  geom_hline(yintercept = 0, alpha = 0.8) +
  geom_point(alpha = 0.8, size = 4, mapping = aes(color = Contrib2)) +
  geom_text_repel(size = 3, point.padding = 0.1) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims2) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  scale_colour_distiller(palette = "YlOrRd", direction = 1) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Main Category Contributors to Axis 2") +
  guides(color = guide_colourbar(title = "Degree of Contribution"))

#Axis 3
plot2.1.3 <- ggplot(tidyvar3,
                mapping = aes(x = Axis2, y = Axis3, label = varnames)) +
  geom_vline(xintercept = 0, alpha = 0.8) +
  geom_hline(yintercept = 0, alpha = 0.8) +
  geom_point(alpha = 0.8, size = 4, mapping = aes(color = Contrib3)) +
  geom_text_repel(size = 3, point.padding = 0.1) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  scale_colour_distiller(palette = "YlOrRd", direction = 1) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Main Category Contributors to Axis 3") +
  guides(color = guide_colourbar(title = "Degree of Contribution"))

#Individuals with contributions above the mean for each principal axis
#Axis 1
plot2.2.1 <- ggplot(tidyind1,
                mapping = aes(x = Axis1, y = Axis2, label = indnames)) +
  geom_vline(xintercept = 0, alpha = 0.8) +
  geom_hline(yintercept = 0, alpha = 0.8) +
  geom_point(alpha = 0.8, size = 4, mapping = aes(color = Contrib1)) +
  geom_text_repel(size = 3, point.padding = 0.1) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims2) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  scale_colour_distiller(palette = "YlOrRd", direction = 1) +
  labs(x = "Dimension 1 (10.23%)",
       y = "Dimension 2 (9.80%)",
       title = "Main Individual Contributors to Axis 1") +
  guides(color = guide_colourbar(title = "Degree of Contribution"))

#Axis 2
plot2.2.2 <- ggplot(tidyind2,
                mapping = aes(x = Axis2, y = Axis3, label = indnames)) +
  geom_vline(xintercept = 0, alpha = 0.8) +
  geom_hline(yintercept = 0, alpha = 0.8) +
  geom_point(alpha = 0.8, size = 4, mapping = aes(color = Contrib2)) +
  geom_text_repel(size = 3, point.padding = 0.1) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  scale_colour_distiller(palette = "YlOrRd", direction = 1) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Main Individual Contributors to Axis 2") +
  guides(color = guide_colourbar(title = "Degree of Contribution"))

#Axis 3
plot2.2.3 <- ggplot(tidyind3,
                mapping = aes(x = Axis2, y = Axis3, label = indnames)) +
  geom_vline(xintercept = 0, alpha = 0.8) +
  geom_hline(yintercept = 0, alpha = 0.8) +
  geom_point(alpha = 0.8, size = 4, mapping = aes(color = Contrib3)) +
  geom_text_repel(size = 3, point.padding = 0.1) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  scale_colour_distiller(palette = "YlOrRd", direction = 1) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Main Individual Contributors to Axis 3") +
  guides(color = guide_colourbar(title = "Degree of Contribution"))

#Plot the supplemental category Twitter influence score against the cloud of categories on axes 1-2
plot3.1.12 <- ggplot() +
  geom_vline(xintercept = 0, color = "gray80") +
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(data = tidyvartable, 
             mapping = aes(x = Axis1, y = Axis2), 
             alpha = 0.5) +
  geom_point(data = twittertable, 
             mapping = aes(x = Axis1, y = Axis2), 
             alpha = 1.0, size = 3, color = "red", shape = 17) +
  geom_text_repel(data = wholetwittertable, 
                  mapping = aes(x = Axis1, y = Axis2, label = varnames), size = 2) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims2) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  labs(x = "Dimension 1 (10.23%)",
       y = "Dimension 2 (9.80%)",
       title = "Supplementary Variables on Principal Axes 1-2",
       subtitle = "Twitter Influence Score")

#Plot the supplemental category Twitter influence score against the cloud of categories on axes 2-3
plot3.1.23 <- ggplot() +
  geom_vline(xintercept = 0, color = "gray80") +
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(data = tidyvartable, 
             mapping = aes(x = Axis2, y = Axis3), 
             alpha = 0.5) +
  geom_point(data = twittertable, 
             mapping = aes(x = Axis2, y = Axis3), 
             alpha = 1.0, size = 3, color = "red", shape = 17) +
  geom_text_repel(data = wholetwittertable, 
                  mapping = aes(x = Axis2, y = Axis3, label = varnames), 
                  size = 2) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Supplementary Variables on Principal Axes 2-3",
       subtitle = "Twitter Influence Score")

#Plot the supplemental category Administration speeches against the cloud of categories on axes 1-2
plot3.2.12 <- ggplot() +
  geom_vline(xintercept = 0, color = "gray80") +
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(data = tidyvartable, 
             mapping = aes(x = Axis1, y = Axis2), 
             alpha = 0.5) +
  geom_point(data = speechtable, 
             mapping = aes(x = Axis1, 
                           y = Axis2, 
                           color = Administration, 
                           shape = Administration), 
             alpha = 1.0, size = 3) +
  geom_text_repel(data = wholespeechtable, 
                  mapping = aes(x = Axis1, y = Axis2, label = varnames), size = 2) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims2) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  labs(x = "Dimension 1 (10.23%)",
       y = "Dimension 2 (9.80%)",
       title = "Supplementary Variables on Principal Axes 1-2",
       subtitle = "Obama vs. Trump Administration Speeches")

#Plot the supplemental category Speeches against the cloud of categories on axes 2-3
plot3.2.23 <- ggplot() +
  geom_vline(xintercept = 0, color = "gray80") +
  geom_hline(yintercept = 0, color = "gray80") +
  geom_point(data = tidyvartable, 
             mapping = aes(x = Axis2, y = Axis3), 
             alpha = 0.5) +
  geom_point(data = speechtable, 
             mapping = aes(x = Axis2, 
                           y = Axis3, 
                           color = Administration, 
                           shape = Administration), 
             alpha = 1.0, size = 3) +
  geom_text_repel(data = wholespeechtable, 
                  mapping = aes(x = Axis2, y = Axis3, label = varnames), size = 2) +
  coord_fixed(ratio = 1) +
  scale_x_continuous(breaks = ticks, limits = lims) +
  scale_y_continuous(breaks = ticks, limits = lims) +
  labs(x = "Dimension 2 (9.80%)",
       y = "Dimension 3 (7.44%)",
       title = "Supplementary Variables on Principal Axes 2-3",
       subtitle = "Obama vs. Trump Administration Speeches")
```